Our initial species selection includes 11 single species groups from the NOBA model. We looked at diet overlap here. Below we will double check survey census outputs against the true values, and then test some survey specifications that may be reasonable for multispecies model testing. Finally we will apply survey functions to the detailed diet data.
fgs <- load_fgs(here("atlantisoutput", "NOBA_March_2020"), "nordic_groups_v04.csv")
lname <- data.frame(Latin = c("*Hippoglossoides platessoides*",
"*Reinhardtius hippoglossoides*",
"*Scomber scombrus*",
"*Melongrammus aeglefinus*",
"*Pollachius virens*",
"*Sebastes mentella*",
"*Micromesistius poutassou*",
"*Clupea harengus*",
"*Gadus morhua*",
"*Boreogadus saida*",
"*Mallotus villosus*"),
Code = c("LRD", "GRH", "MAC", "HAD", "SAI", "RED",
"BWH", "SSH", "NCO", "PCO", "CAP")
)
sppsubset <- merge(fgs, lname, all.y = TRUE)
spptable <- sppsubset %>%
arrange(Index) %>%
select(Name, Long.Name, Latin)
knitr::kable(spptable, col.names = c("Model name", "Full name", "Latin name"))
| Model name | Full name | Latin name |
|---|---|---|
| Long_rough_dab | Long rough dab | Hippoglossoides platessoides |
| Green_halibut | Greenland halibut | Reinhardtius hippoglossoides |
| Mackerel | Mackerel | Scomber scombrus |
| Haddock | Haddock | Melongrammus aeglefinus |
| Saithe | Saithe | Pollachius virens |
| Redfish | Redfish | Sebastes mentella |
| Blue_whiting | Blue whiting | Micromesistius poutassou |
| Norwegian_ssh | Norwegian spring spawning herring | Clupea harengus |
| North_atl_cod | Northeast Atlantic cod | Gadus morhua |
| Polar_cod | Polar cod | Boreogadus saida |
| Capelin | Capelin | Mallotus villosus |
Note that survey config files must now name the survey.
NOBA2config.R looks like this (adjusted from Alfonso’s original):
d.name <- here("atlantisoutput","NOBA_march_2020")
functional.groups.file <- "nordic_groups_v04.csv"
biomass.pools.file <- "nordic_biol_v23.nc"
biol.prm.file <- "nordic_biol_incl_harv_v_007_3.prm"
box.file <- "Nordic02.bgm"
initial.conditions.file <- "nordic_biol_v23.nc"
run.prm.file <- "nordic_run_v01.xml"
scenario.name <- "nordic_runresults_01"
bioind.file <- "nordic_runresults_01BiomIndx.txt"
catch.file <- "nordic_runresults_01Catch.txt"
annage <- TRUE
fisheries.file <- "NoBAFisheries.csv"
omdimensions.R standardizes timesteps, etc. (this is part of atlantisom and should not need to be changed by the user):
#survey species inherited from omlist_ss
survspp <- omlist_ss$species_ss
# survey season and other time dimensioning parameters
# generalized timesteps all models
noutsteps <- omlist_ss$runpar$tstop/omlist_ss$runpar$outputstep
timeall <- c(0:noutsteps)
stepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutinc
midptyr <- round(median(seq(0,stepperyr)))
# model areas, subset in surveyconfig
allboxes <- c(0:(omlist_ss$boxpars$nbox - 1))
# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc
# survey selectivity (agecl based)
sp_age <- omlist_ss$funct.group_ss[, c("Name", "NumCohorts", "NumAgeClassSize")]
# should return all age classes fully sampled (Atlantis output is 10 age groups per spp)
n_age_classes <- sp_age$NumCohorts
# changed below for multiple species NOTE survspp alphabetical; NOT in order of fgs!!
# this gives correct names
age_classes <- sapply(n_age_classes, seq)
names(age_classes)<-sp_age$Name
n_annages <- sp_age$NumCohorts * sp_age$NumAgeClassSize
# changed below for multiple species
annages <- sapply(n_annages, seq)
names(annages)<-sp_age$Name
mssurvey_spring.R and mssurvey_fall.R configure the fishery independent surveys (in this census test, surveys sample all model polygons in all years and have efficiency of 1 for all species, with no size selectivity):
# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity
# Survey name
survey.name="BTS_spring_allbox_effic1"
#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5
#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May,
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)
survey_sample_time <- 1 # spring survey
#The last timestep to sample
total_sample <- noutsteps-1 #495
#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
total_sample, by=timestep)
survtime <- survey_sample_full
# survey area
# should return all model areas
survboxes <- allboxes
# survey efficiency (q)
# should return a perfectly efficient survey
surveffic <- data.frame(species=survspp,
efficiency=rep(1.0,length(survspp)))
# survey selectivity (agecl based)
# this is by age class, need to change to use with ANNAGEBIO output
#survselex <- data.frame(species=rep(names(age_classes), each=n_age_classes),
# agecl=rep(c(1:n_age_classes),length(survspp)),
# selex=rep(1.0,length(survspp)*n_age_classes))
# for annage output uses names(annages) NOT alphabetical survspp
survselex <- data.frame(species=rep(names(annages), n_annages), #
agecl=unlist(sapply(n_annages,seq)),
selex=rep(1.0,sum(n_annages)))
# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
surveffN <- data.frame(species=survspp, effN=rep(100000, length(survspp)))
# survey index cv needed for sample_survey_xxx
# cv = 0.1
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))
# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1
# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200
# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity
# Survey name
survey.name="BTS_fall_allbox_effic1"
#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5
#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May,
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)
survey_sample_time <- 3 # fall survey
#The last timestep to sample
total_sample <- noutsteps-1 #495
#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
total_sample, by=timestep)
survtime <- survey_sample_full
# survey area
# should return all model areas
survboxes <- allboxes
# survey efficiency (q)
# should return a perfectly efficient survey
surveffic <- data.frame(species=survspp,
efficiency=rep(1.0,length(survspp)))
# survey selectivity (agecl based)
# this is by age class, need to change to use with ANNAGEBIO output
#survselex <- data.frame(species=rep(survspp, each=n_age_classes),
# agecl=rep(c(1:n_age_classes),length(survspp)),
# selex=rep(1.0,length(survspp)*n_age_classes))
# for annage output
survselex <- data.frame(species=rep(names(annages), n_annages), #
agecl=unlist(sapply(n_annages,seq)),
selex=rep(1.0,sum(n_annages)))
# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
surveffN <- data.frame(species=survspp, effN=rep(100000, length(survspp)))
# survey index cv needed for sample_survey_xxx
# cv = 0.1
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))
# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1
# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200
msfishery.R configures the fishery dependent data:
# Default fishery configuration here is a census
# Inherits species from input omlist_ss
fishspp <- omlist_ss$species_ss
#Number of years of data to pull
nyears <- 50
#Atlantis initialization period in years
burnin <- 30
# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc
# same time dimensioning parameters as in surveycensus.R
#Vector of indices of catch in numbers to pull (by timestep to sum)
fish_sample_full <- c(0:total_sample) #total_sample defined in sardinesurvey.R
fish_burnin <- burnin*fstepperyr+1
fish_nyears <- nyears*fstepperyr
fish_times <- fish_sample_full[fish_burnin:(fish_burnin+fish_nyears-1)]
fish_timesteps <- seq(fish_times[fstepperyr], max(fish_times), by=fstepperyr) #last timestep
#fish_years <- unique(floor(fish_times/fstepperyr)+1) # my original
fish_years <- unique(floor(fish_times/fstepperyr)) #from Christine's new sardine_config.R
fishtime <- fish_times
# fishery sampling area
# should return all model areas, this assumes you see everything that it caught
fishboxes <- c(0:(omlist_ss$boxpars$nbox - 1))
# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
fisheffN <- data.frame(species=survspp, effN=rep(1000, length(survspp)))
#this adjusts for subannual fishery output so original effN is for whole year
fisheffN$effN <- fisheffN$effN/fstepperyr
# fishery catch cv can be used in sample_survey_biomass
# perfect observation
fish_cv <- data.frame(species=survspp, cv=rep(0.01,length(survspp)))
These survey and fishery configurations were used previously to generate saved census survey output.
Using read_savedsurvs() we get back the same objects generated by the wrapper functions om_index() and om_comps().
source(here("config", "NOBA2config.R"))
omlist_ss <- readRDS(file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))
source(here("config/omdimensions.R"))
survObsBiom <- read_savedsurvs(d.name, 'survB')
age_comp_data <- read_savedsurvs(d.name, 'survAge')
len_comp_data <- read_savedsurvs(d.name, 'survLen')
wtage <- read_savedsurvs(d.name, 'survWtage')
annage_comp_data <- read_savedsurvs(d.name, 'survAnnAge')
annage_wtage <- read_savedsurvs(d.name, 'survAnnWtage')
# no function needed for saved fishery output as long as there is only one fishery specification file
catchbio_ss <- readRDS(file.path(d.name, paste0(scenario.name, "fishCatch.rds")))
catchbio_ss <- catchbio_ss[[1]] #in same time unit as surveys
fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
fish_age_comp <- fish_age_comp[[1]]
#add this to om_indices function so that this has years when read in
#fish_age_comp$time <- fish_age_comp$time/fstepperyr
fish_len_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsLenComp.rds")))
#len_comp_data <- len_comp_data[[1]]
fish_len_comp_data <- fish_len_comp_data[[1]]
#add this to om_indices function so that this has years when read in
#fish_len_comp_data$time <- as.integer(floor(fish_len_comp_data$time/fstepperyr))
fishwtage <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsWtAtAge.rds")))
fishwtage <- fishwtage[[1]] #need to change time units?
fish_annage_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsFullAgeComp.rds")))
fish_annage_comp <- fish_annage_comp[[1]]
#add this to om_indices function so that this has years when read in
#fish_annage_comp$time <- fish_annage_comp$time/fstepperyr
fish_annage_wtage <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsFullWtAtAge.rds")))
fish_annage_wtage <- fish_annage_wtage[[1]] #need to change time units?
Crunch a few things to get true age comps
# age classes
totNagecl <- omlist_ss$truenums_ss %>%
group_by(species, agecl, time) %>%
summarize(numAtAge = sum(atoutput))
totN <- totNagecl %>%
group_by(species, time) %>%
summarize(totN = sum(numAtAge))
trueNagecl <- merge(totNagecl, totN)
# annual ages
totNage <- omlist_ss$truenumsage_ss %>%
group_by(species, agecl, time) %>%
summarize(numAtAge = sum(atoutput))
#this should be exactly the same as totN but just making sure...
totN2 <- totNage %>%
group_by(species, time) %>%
summarize(totN = sum(numAtAge))
trueNage <- merge(totNage, totN2)
# fishery catch at age class
totCagecl <- omlist_ss$truecatchnum_ss %>%
group_by(species, agecl, time) %>%
summarize(numAtAge = sum(atoutput))
totCAAcl <- totCagecl %>%
group_by(species, time) %>%
summarize(totN = sum(numAtAge))
trueCAAcl <- merge(totCagecl, totCAAcl)
# fishery catch at annual age
totCage <- omlist_ss$truecatchage_ss %>%
group_by(species, agecl, time) %>%
summarize(numAtAge = sum(atoutput))
totCAA <- totCage %>%
group_by(species, time) %>%
summarize(totN = sum(numAtAge))
trueCAA <- merge(totCage, totCAA)
My plotting functions: too specialized to include in the package?
# plot biomass time series facet wrapped by species
plotB <- function(dat, truedat=NULL){
ggplot() +
geom_line(data=dat, aes(x=time/stepperyr,y=atoutput, color="Survey Biomass"),
alpha = 10/10) +
{if(!is.null(truedat)) geom_line(data=truedat, aes(x=time/365,y=atoutput, color="True B"), alpha = 3/10)} +
theme_tufte() +
theme(legend.position = "top") +
xlab("model year") +
ylab("tons") +
labs(colour=scenario.name) +
facet_wrap(~species, scales="free")
}
# make a catch series function that can be split by fleet? this doesnt
# also note different time (days) from model timestep in all other output
plotC <- function(dat, truedat=NULL){
ggplot() +
geom_line(data=dat, aes(x=time/365,y=atoutput, color="Catch biomass"),
alpha = 10/10) +
{if(!is.null(truedat)) geom_line(data=truedat, aes(x=time/365,y=atoutput, color="True Catch"), alpha = 3/10)} +
theme_tufte() +
theme(legend.position = "top") +
xlab("model year") +
ylab("tons") +
labs(colour=scenario.name) +
facet_wrap(~species, scales="free")
}
# plot length frequencies by timestep (one species)
plotlen <- function(dat, truedat=NULL){
ggplot(dat, aes(upper.bins)) +
geom_bar(aes(weight = atoutput)) +
theme_tufte() +
xlab("length (cm)") +
ylab("number") +
labs(subtitle = paste(scenario.name,
dat$species)) +
facet_wrap(~time, ncol=6, scales="free_y")
}
# plot numbers at age by timestep (one species)
Natageplot <- function(dat, effN=1, truedat=NULL){
ggplot() +
geom_point(data=dat, aes(x=agecl, y=atoutput/effN, color="Est Comp")) +
{if(!is.null(truedat)) geom_line(data=dat, aes(x=agecl, y=numAtAge/totN, color="True Comp"))} +
theme_tufte() +
theme(legend.position = "bottom") +
xlab("age/agecl") +
ylab("number") +
labs(subtitle = paste(scenario.name,
dat$species)) +
facet_wrap(~time, ncol=6, scales="free_y")
}
# plot weight at age time series facet wrapped by species
wageplot <- function(dat, truedat=NULL){
ggplot(dat, aes(time/stepperyr, atoutput)) +
geom_line(aes(colour = factor(agecl))) +
theme_tufte() +
theme(legend.position = "bottom") +
xlab("model year") +
ylab("average individual weight (g)") +
labs(subtitle = paste0(scenario.name)) +
facet_wrap(c("species"), scales="free_y")
}
We are using the full time series for biomass and weight at age, and a subsample of years 30-53 for compositions. In this test we see the effect of surveying in one season, with migratory species missed or only partially represented in some surveys, and with the overall surveyed biomass being at the lower or higher end of the range depending on the time of year. (Observation error is included in the survey biomass time series.) Samples for length and age assume a total of 100000 fish were randomly sampled for each species in each survey (surveffN defined in the survey config files).
# compare with true output (all timesteps)
for(s in names(survObsBiom)){
cat(" \n##### ", s," \n")
print(plotB(survObsBiom[[s]][[1]], omlist_ss$truetotbio_ss))
cat(" \n")
}
# plots survey only
# for(s in names(survObsBiom)){
# print(plotB(survObsBiom[[s]][[1]]))
# }
for(s in names(len_comp_data)){
cat(" \n##### ", s," \n")
lcompsub <- as.data.frame(len_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
group_by(species) %>%
group_map(~ plotlen(.x), keep = TRUE)
for(i in 1:length(lcompsub)) {
print(lcompsub[[i]])
}
cat(" \n")
}
for(s in names(age_comp_data)){
cat(" \n##### ", s," \n")
acompsub <- as.data.frame(age_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
group_by(species) %>%
left_join(., trueNagecl) %>%
#group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
group_map(~ Natageplot(.x, effN = 100000, truedat = 1), keep = TRUE) # plots merged true age comp
for(i in 1:length(acompsub)) {
print(acompsub[[i]])
}
cat(" \n")
}
for(s in names(annage_comp_data)){
cat(" \n##### ", s," \n")
acompsub <- as.data.frame(annage_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
group_by(species) %>%
left_join(., trueNage) %>%
#group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
group_map(~ Natageplot(.x, effN = 100000, truedat = 1), keep = TRUE) # plots merged true age comp
for(i in 1:length(acompsub)) {
print(acompsub[[i]])
}
cat(" \n")
}
for(s in names(wtage)){
cat(" \n##### ", s," \n")
print(wageplot(wtage[[s]][[1]]))
cat(" \n")
}
# annage_wtage read in above
# annage_wtage <- annage_wtage[[1]] #this still has second list component, diagnostic plot
# annage_wtage <- annage_wtage[[1]] #still in this old version but just removed from function.
for(s in names(annage_wtage)){
cat(" \n##### ", s," \n")
print(wageplot(annage_wtage[[s]][[1]][[1]]))
cat(" \n")
}
We are using the full time series for total catch and weight at age, and a subsample of years 30-35 for compositions. Fishery length and age composition is sampled every output timestep. Fishery effN for length and age sampling was 1000 fish per year (200 per time step), so this is reflected in the composition plots.
# compare with true catch, should match
plotC(catchbio_ss, omlist_ss$truecatchbio_ss)
lcompsub <- as.data.frame(fish_len_comp_data) %>% filter(time %in% c(150:175)) %>%
group_by(species) %>%
group_map(~ plotlen(.x), keep = TRUE)
for(i in 1:length(lcompsub)) {
print(lcompsub[[i]])
}
acompsub <- as.data.frame(fish_age_comp) %>% filter(time %in% c(150:175)) %>%
group_by(species) %>%
left_join(., trueCAAcl) %>%
#group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
group_map(~ Natageplot(.x, effN = 200, truedat = 1), keep = TRUE) # plots merged true age comp
for(i in 1:length(acompsub)) {
print(acompsub[[i]])
}
acompsub <- as.data.frame(fish_annage_comp) %>% filter(time %in% c(150:175)) %>%
group_by(species) %>%
left_join(., trueCAA) %>%
#group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
group_map(~ Natageplot(.x, effN = 200, truedat = 1), keep = TRUE) # plots merged true age comp
for(i in 1:length(acompsub)) {
print(acompsub[[i]])
}
wageplot(fishwtage)
#diagnostic plot component still in this old version but just removed from function.
wageplot(fish_annage_wtage[[1]])